Summarise taxonmic distribution in the EMP dataset

ftp <- "ftp://ftp.microbio.me/emp/release1/"
biom_file <- "emp_cr_silva_16S_123.subset_10k.biom"
map_file <- "emp_qiime_mapping_subset_10k.tsv"
new_map_file <- "emp_qiime_mapping_subset_10k_noBLNK.csv"
Rank <- "Phylum"
Grouping_var <- "empo_3"
EMPO_levels2keep <- c("Soil (non-saline)", "Plant rhizosphere", "Animal corpus", "Water (saline)", "Sediment (saline)", "Water (non-saline)", "Hypersaline (saline)", "Sediment (non-saline)", "Animal distal gut", "Animal proximal gut", "Aerosol (non-saline)", "Animal secretion", "Animal surface", "Plant corpus", "Surface (non-saline)", "Plant surface", "Surface (saline)") # remove unwanted environments
s_size = 1000 # genomes to sample per community

biom_ftp <- paste0(ftp, "otu_tables/closed_ref_silva/", biom_file)
map_ftp <- paste0(ftp, "mapping_files/", map_file)
download.file(biom_ftp, biom_file)
download.file(map_ftp, map_file)

# fill-in missing values with NA, convert to CSV
system("awk 'BEGIN { FS = OFS = \"\t\" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = \"NA\" }; 1' emp_qiime_mapping_subset_10k.tsv  |  awk 'BEGIN { FS=\",\"; OFS=\";\" } {$1=$1; print}' |  awk 'BEGIN { FS=\"\t\"; OFS=\",\" } {$1=$1; print}' > emp_qiime_mapping_subset_10k_noBLNK.csv")

(EMP_otu_tax <- import_biom(biom_file, parseFunction = parse_taxonomy_silva_128, 
                            parallel = TRUE))

read_csv(new_map_file,
                    trim_ws = TRUE) %>% 
  column_to_rownames("#SampleID") %>% 
  sample_data() ->
  EMP_map

EMP <- merge_phyloseq(EMP_otu_tax, EMP_map)

# for debugging and testing
# EMP %>%
#   subset_samples(empo_3 %in% EMPO_levels2keep) %>%
#   prune_samples(sample_names(.)[sample(1:nsamples(.), 100)], .) %>%
#   phyloseq::filter_taxa(., function(x) sum(x) > 0, TRUE) %>%
#   prune_taxa(taxa_names(.)[sample(1:ntaxa(.), 1000)], .) ->
#   EMP_subset

EMP %>% 
  PS_merge_samples(grouping_name = Grouping_var, fun = "mean") ->
  EMP_merge

save(EMP, file = paste0("EMP.RDS"))
save(EMP_merge, file = paste0("EMP_merge.RDS"))
# load("./EMP.RDS")
# load("./EMP_merge.RDS")

# summarise abundance by environment and phylum
EMP %>%
  tax_glom(taxrank = Rank) %>%                     # agglomerate at 'Rank' level
  psmelt() %>%                                         # Melt to long format
  arrange(Rank) %>%                                  # arrange by 'Rank'
  group_by_at(c(Grouping_var, Rank)) %>% 
  # summarise and create a column with the relative abundance
  summarise(`Abundance (mean)` = mean(Abundance),
            n = n(),  # calculates the sample size per group
            SD = sd(Abundance, na.rm = T)) %>%  # calculates the standard error of each group
  mutate(Freq = (100 * `Abundance (mean)` / sum(`Abundance (mean)`)),
         `Freq SD` = (100 * SD / sum(`Abundance (mean)`))) -> # abundance per empo_3 cat
  taxa_sum_df

# determine the phylum order by decreasing abundance (for plotting)
EMP_merge %>%
  tax_glom(taxrank = Rank) %>%                     # agglomerate at 'Rank' level
    transform_sample_counts(., function(x) x / sum(x) * 100) %>% 
  psmelt() %>%                                        # Melt to long format
  arrange(Rank)  %>%                                  # arrange by 'Rank'
  group_by_at(c("OTU", Rank)) %>% 
  summarise(`Abundance (mean)` = mean(Abundance)) %>%
  arrange(`Abundance (mean)`) ->
  Rank_order

# print summary table
taxa_sum_df %>% 
  filter(`Abundance (mean)` != 0) %>%
  kable(., digits = c( 0, 1)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)